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We formulate a method to study two-body correlations in a condensate of identical bosons. 
We use the adiabatic hyperspheric approach and assume a Faddeev like decomposition of the wave 
function. We derive for a fixed hyperradius an integro-differential equation for the angular eigenvalue 
and wave function. We discuss properties of the solutions and illustrate with numerical results. The 
interaction energy is for A'^ ~ 20 five times smaller than that of the Gross-Pitaevskii equation. 
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Introduction. Few-body correlations often express the 
distinguishing characteristic features of an iV-body sys- 
tem Q . Two-body correlations are not only the simplest 
but in most cases also the most important. Higher order 
correlations have a tendency to either strongly confine 
spatially or correlate into clusters of particles effectively 
reducing the correlations to lower order. Exceptions are 
the three-body correlations decisive for the stability of 
Borromcan systems known for dripline nuclei Q . 

Non-correlated mean-field computations of nuclei with 
the free two-body nucleon-nucleon interaction produce 
disastrously wrong results. Two-body correlations com- 
pensating for the short range hard core repulsion are ab- 
solutely necessary [|j . For atoms the effective interaction 
is also strongly repulsive at shorter distances. Further- 
more for an atomic Bose condensate a short-range two- 
body attraction produces diatomic recombination and 
thereby the atoms decay out of the condensate For 
both nuclei and molecules correlations are decisive. For 
nuclei various methods have been designed to deal with 
this problem, i.e. Jastrow theory, Bruckner theory, effec- 
tive interactions in model spaces 

For Bose condensates the Gross-Pitaevskii (non- 
correlated) mean-field equation has been the starting 
point since the first observation of the condensate in 1995 
1^ . The wave function does not include any correlations 
and the assumed repulsive (5-interaction has the imme- 
diate consequence that the short range behavior can- 
not be described even qualitatively correct. Repulsive 
zero range potentials are not physically meaningful |^] 
although sometimes useful in selected Hilbert spaces. At- 
tractive (5-interactions in three dimensions lead to diver- 
gencies demanding renormalization or a change of 
boundary conditions [|o|. When two-body bound states 
appear even dimer condensates may be possible ]ll| ] . 

To include correlations we must necessarily go beyond 
the mean-field Hartree-Fock-Boguliubov approximation. 
Then finite range potentials with realistic features can 
as well be used as the starting point in the theoretical 
formulation. The other crucial ingredient is the degrees 
of freedom or, equivalently, the Hilbert space. An in- 
teresting formulation was recently introduced in terms 
of generalized hyperspherical coordinates and an adia- 
batic expansion with the hyperradius as the adiabatic 



coordinate [Q. Still only a zero-range interaction was 
used with the lowest hyperspherical angular wave func- 
tion. These are very crude approximations, since the ex- 
pansion in the hyperspherical basis necessarily must be 
extremely slowly converging for large(r) hyperradii. 

The purpose of this letter is to go a step further and 
establish the optimum equation to determine two-body 
correlations. We shall derive a practical and realistic 
equation applicable to a system of N identical bosons. 
We shall allow general two-body interactions and in par- 
ticular both attractive and repulsive finite range poten- 
tials. The resulting equation may be considered as an 
alternative to the Gross-Pitaevskii equation, but allow- 
ing two-body correlations. 

Theory. The system of N identical interacting bosons 
of mass m may be described by the coordinates rj or CM 
(R) and the A; = 1,2,...,7V— 1 Jacobi coordinates 
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The hyperspherical coordinates are given by the hyper- 
radius p and the hyperangles € [0,7r/2] 
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sin aij , and we distinguish be- 



where rij = {ri — fj 
tween one and two indices on a. The remaining 2{N — 1) 
angles are the directions of the — 1 ?7fc-vectors. All the 
angles together, including the a's, are denoted fl. 

Removing the center of mass motion the intrinsic 
hamiltonian H for a trap of angular frequency lo becomes 
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where Vij = V{rij) is the two-body interaction and the 



internal kinetic energy T is of the form |13 



T 



A^ = - 



2m 



„3N- 



^p3N-4_^ 3_^2 

dp dp p^ 



3Ar-9- (3A-5)cos2a d 
sin 2a da 



D 



(4) 
(5) 



1 



where D only contains derivatives with respect to angles 
different from a = an- The hamiltonian is then 
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where Vij — Vij2mp^ /h is a dimensionless potential. 

The total wave function ^ , obeying H"^ = E'if, is for 
each p expanded as 
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where fn{p) are the expansion coefhcients on the solu- 
tions $„(p, J7) to the angular eigenvalue equation 



(/in-A„)$„(p,r!) = 0. 
The coupled set of radial equations are then ||l3| 
Kipl , 2m{U{p)-E y 

E(20in'(P)^ + Qw(p))/"'' 
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with U{p) from the external trap and centrifugal barrier 
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as 



2 ' 8m/92 
Then is decomposed in Faddeev components 

JV 

$„(p,i])= 4;\p,^)- (13) 

t<j=i 

We assume that each two-particle amplitude 4>ij can be 
described by s-waves and therefore only depends on p 
and the distance rij |Q. This does not exclude higher 
partial waves in the total wave function, e.g. expressing 
the "12" s-wave component in relative "34" -coordinates 
requires non-zero angular momenta. The boson symme- 
try implies that fi) = (j)'^"-\p,aij) = ^(ay) are 
identical functions of the different coordinates a. 



The angular eigenvalue A in eq.(|9|) is given by 
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We now vary (j)^2 for each p and fi nd the condition for 
extremum of A. From eqs.(^) and (|lj) we then obtain 
the integro-differential equation 
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P + v{a)-X] 0(a) + / dTG{T, a) = 



The integrals in eq.(|15|) involve at most six particle 
coordinates. By appropriate choices [H) of Jacobi sys- 
tems these can be expressed in terms of the five vectors 
77JV-1, . . . , ?7)v-5. One variable (a = au) is the argument 
of the variational function 0i2 and not an integration 
variable. Furthermore 7 variables only enter in the phase 
space and can be integrated analytically leaving at most 
five-dimensional integrals. With v{a) = vi2{V2psma) 
we rewrite eq.(E3) 
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where r denotes the remaining 5 variables and G is a 
definite homogeneous linear combination of (j){aik{T, a)). 

The integrals in eq.(|l6|) reduce to at most two dimen- 
sions when the two-body interaction-range h is much 
less than the characteristic length p of the system. 
To illustrate we use a gaussian potential V{rij) = 
Vq exp(— r|^/62), but the results depend mainly on = 

^/7rm6'^Vo/4?i^, i.e. the Born-approximation of the scat- 
tering length Qs- When mb^Vo <SC we have Up « Og. 
For 6 <C p we get, almost independent of the potential. 
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Here T is the gamma function, sin /3 — tan Q!/\/3, 6 is the 
truth function, and R, a functional of v and cj) and a func- 
tion of a, arises from rotations of interaction and wave 
function from one set of Jacobi coordinates to another. 

The ^-interaction with a constant wave function, la- 
beled K = {), leads to the eigenvalue |12 



which for large N coincides with the term v{a = 0) in 
eq.(|l^). The interaction dependent part of the strength 
is collected in the parameter ap. 

The limit of A^^o is not, even for large N, obtained by 
solving eqs.(|l^) and (|l7|), since the highest power of N 
is found in one of the rotated terms in R proportional to 
N'^v{a). In the zero-range approximation the two non- 
local terms in eq.(pT|) proportional to (j>{Q) are very sen- 
sitive to the initial two-body potential v{a) which is very 
large and of very short range already for moderate values 
of p. These terms are therefore crucial for the develope- 
ment of correlations. 
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Quantitative behavior. We now solve numerically 
eq.(|l^) using eq.(|l7|) with the realistic parameter set in 
| p2t for the condensate of ^^i?&-atoms. For numerical il- 
lustration we choose N — 10, 20, 30. 

The lowest angular eigenvalues for iV = 20 are shown 
in Fig. |l| in the relevant range of hyperradii selected by 
the external trap. The minimum of U{p) is located near 
p ~ bt\/3N/2 where bt = y/h/(jmj) is the trap length. 
Inclusion of the interaction terms, expressed in eq.(|l^) by 
A through the effective radial potential U{p) + h'^ \/2mp^ , 
push the minimum outwards. For the repulsive interac- 
tion the angular eigenvalues decrease with p due to an 
increasing average two-body distance. 

We first only include the local terms (oc 0(q;)) in 
eq.(|l^). The angular energy in Fig. ^ is then larger than 
the Xk=o value, because these terms are repulsive. In- 
cluding also the i?-terms in eq. leads to rather similar 
results, where however, the correlations are exploited giv- 
ing a lower energy. Finally, inclusion of all terms reduces 
the angular energy to only 19% of the Xk=o value. 
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FIG. 1. The lowest angular eigenvalue A obtained from 
eq.(|l^) as function of p for 20 *^i?6-atoms. The parame- 
ters are ap = 5.29 nm, b = Sp/lO, v = uj/{2n) = 200 Hz, 
{bt = 763 nm). The decreasing curves show the energ y i n 
eq.(po|), the energy including only the local terms of eq.(p^, 
local plus i?-terms (rotated) and the full energy with all terms. 
Shown is also U{p) (parabolic curve) in eq.([L3) along with the 
free center of mass energy 3(A'^ — l)?ia;/2 (horizontal line). 

The corresponding angular wave functions in Fig. |^ are 
all zero at a = due to the vanishing volume clement. 
The almost discontinous behavior at this point of the full 
and rotated angular wave function is caused by the very 
short range of the initial repulsive two-body potential 
v{a). The K — wave function has no nodes and the lo- 
cal terms maintain this structure but marginally shifted 
to larger a- values by the repulsive potential. The rotated 
terms introduce one node indicating the substantial re- 
structuring due to correlations. The higher kinetic en- 
ergy is more than compensated by the correlations build 



up to avoid the repulsion at short distance. The full so- 
lution maintains qualitatively this behavior, but now the 
probability is shifted to larger a or equivalently to larger 
distances between each pair of particles. 

The terms proportional to (/)(0) decreases the angular 
energy drastically. This is consistent with the node in the 
wave function which lowers the angular energy because 
then large \<p{0) \ results in smaller A. More oscillations al- 
low constructive interference which in turn substantially 
lowers the potential energy. The result is the increasing 
deviation of A from \k=o for decreasing p, see Fig. ^ 
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FIG. 2. The wave function (p multiplied by the square root 

-,s{3iv-7)/2 function 

/20 • 559 nm) for the A's in Fig. 0. 



of the volume element, sinacos'''^ '^'''^ a, as function of a 



for p = 2500 nm (= 



The N and p dependences of the angular wave func- 
tions are shown in Fig. ||. The ^/N scaling follow the 
N dependence of the potential minimum of U (p) , see 
Fig. ^ The strong variation for small a for iV = 30 dis- 
appears with increasing p. The same tendency, but less 
pronounced, is found for N ^ 10. The size of the an- 
gular wave functions for a ~ decreases with increasing 
p. In the relevant parameter interval defined by the trap 
the p-dependence of the angular wave function is weak 
except at small a corresponding to distances inside the 
two-body potential. The corresponding angular eigen- 
value A relative to Ai^^o is essentially independent of p, 
i.e. 0.33, 0.19 and 0.16 for N = 10,20 and 30, respec- 
tively. This confirms numerically our conjecture that the 
K — behavior obtained in with a (5-interaction is 
not approached with increasing N at large distances. 

The total energies obtained by solving the radial equa- 
tion in eq.(p^) reflect the relative sizes of the angular 
eigenvalues in Fig. ^ The radial wave functions fn{p) 
resemble the K = solutions. However, it is essential 
to understand that even if the effective radial potential 
only deviates insignificantly from the K — potential the 
corresponding angular wave functions may differ enor- 
mously, i.e. none versus many oscillations. Thus effects 
of correlations may easily be strong without showing up 
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in the total energy and perhaps also not even significantly 
in the interaction energy. 

The Gross-Pitaevskii mean-field solution does not re- 
move the (small) spurious contribution from the center of 
mass motion. Still we compare directly and find that our 
total energy only is slightly smaller than this total mean- 
field energy. However, our interaction energy, exclud- 
ing the external trap energy, is about 3, 5, and 7 times 
smaller for N — 10, 20, and 30, respectively. This reflects 
the huge difference in structure between the mean-field 
and correlated solutions. 
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FIG. 3. Tlie same as Fig. ^ witli full solutions for hyper- 
radii p = p/ViV = 100 nm, 400 nm, 6400 nm and iV = 10, 30. 



two-body correlations. In the mean-field approximation 
these correlations are by definition completely absent and 
the average distance between pairs of particles is there- 
fore very small. The accurate detailed behavior can only 
be obtained by explicit inclusion of correlations. 

Perspectives. In conclusion, we have variationally de- 
rived a linear one-dimensional intcgro-differential equa- 
tion describing N identical bosons in a trap. The equa- 
tion involves five-dimensional integrals, a tremendous 
simplification from the original A^-body problem. We 
assume the Faddeev angular decomposition of the wave 
function and use the hyperradius as the adiabatic coor- 
dinate. A further reduction to at most two-dimensional 
integrals is achieved for short range interactions for the 
non-local parts of the equation. The interaction may 
be attractive and of finite range with bound states and 
with both signs of the scattering length. This opens the 
possibility of realistic computation of the diatomic re- 
combination rate which in the present terminology is a 
process originating from the first excited adiabatic state 
(condensate) ending in the lowest adiabatic (diatomic 
bound) state. This process is essential for the stability of 
atomic Bose condensates. Our equation is an alternative 
to the Gross-Pitaevskii equation, but apparently both 
features of mean-field and few-body correlations are now 
included in a unified approach. Other interesting per- 
spectives are extensions to three-body correlations, to 
one and two space dimensions, to non-identical bosons, 
to fcrmion systems and to systems of mixed symmetry. 
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particles is = 
nean square "radius" 



((/)(a)|sin a|(/)(a)) obtained with normalized (f), quanti- 
fies this distance for each p. For the p-independent K = 
wave function we find 1/{N — 1) analytically whereas our 
full wave function with the correlations for = 20 gives 
the much larger distance of about 2A/{N — 1) roughly 
independent of N and p, see Figs. I and |. The full mean 
square radius is now obtained by multipliplying (sin^ a) 
with 2(p^) K, 3{N — 1)6( found as the minimum of U in 
Fig. |l] or as the harmonic oscillator expectation value. 
Thus by increasing N the decreasing two-particle aver- 
age distance for fixed p is precisely compensated by the 
increasing position of the minimum of the potential. 

The result is (r^./fo^) = 3 for X = and 6.83, 7.30 
and 7.31 for the full wave function for N = 10, 20, 30, re- 
spectively. These values can be compared to the Gross- 
Pitaevskii results through the harmonic oscillator rela- 
tion {N - 1)(4) = 2iV((rf) - (i?^)), where iV(i?2) ^ 
36j/2 is the center of mass expectation value for the 
mean square radius. With the mean-field result numeri- 
cally computed we obtain {rfj/bf) « 3.1 almost indepen- 
dent of N, i.e. almost equal to the K = result but 
much smaller than for the correlated solution. Thus the 
K ~ assumption do not even approximately describe 
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